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We study the interplay between the fluid-crystal transition and the glass transition of elastic sphere 
system with polydispersity using nonequilibrium molecular dynamics simulations. It is found that 
the end point of the crystal-fluid transition line, which corresponds to the critical polydispersity 
above which the crystal state is unstable, is on the glass transition line. This means that crystal 
and fluid states at the melting point becomes less distinguishable as polydispersity increases and 
finally they become identical state, i.e., marginal glass state, at critical polydispersity. 



I. INTRODUCTION 

Recently, jammed (amorphous) solid has attracted 
more attention in two aspects; 1) the origin of rigidity 
which fluid lacks, and 2) whether some essential differ- 
ences exists between the jammed solid and crystal except 
for the positional order. The former corresponds to the 
long standing problem of glass transition |l[ and jam- 
ming transition m, and the latter is related to the dense 
packing problem |3[ . Since whether the solid has or lacks 
the positional order is complexly related to properties of 
materials, it is difficult to find which factor is important. 
Systematic study is, however, possible in polydisperse 
particle system, which consists of multiple ingredients 
with various sizes or shapes. It is empirically known that 
binary fluid mixture, with more than 10% size-dispersion, 
often exhibits glass transition Q while monodisperse sim- 
ple fluid exhibits crystallization transition. Therefore, 
polydispersity is one of the most important factors for 
the glass transition. While the system may involve the 
glass transition for high enough polydispersity and the 
simple crystallization for low enough polydispersity, there 
remains a gap in knowledge between these two regimes. 
In order to study the effect of polydispersity, accurate 
control of polydispersity is required, which is difficult in 
experiments. 

Here, we consider the polydisperse hard-sphere (HS) 
system, which is one of the simplest models to exhibit 
both of fluid and crystal phases. While it is well known 
that monodisperse HS system takes a first order melt- 
ing/crystallization transition, so-called the Alder transi- 
tion, by increasing/decreasing pressure or density [1,0], 
much attention recently has been attracted to the prob- 
lem; how polydispersity affects this transition. Most re- 
markable finding is that the discontinuity at the melting 
point, i.e., density gap between the fluid and the crys- 
tal, decreases as the strength of polydispersity increases 
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FIG. 1: (color online) Phase diagram showing the polydisper- 
sity vs packing fraction (or pressure, inset) plane. There are 
three equilibrium phases: fluid, crystal, fluid-crystal coexis- 
tence. The boundary between the fluid and the disordered 
solid indicates a dynamical transition. 



and finally vanishes at a certain critical point [7H2( (see 
Fig- HI- This is similar to the well-known liquid-gas crit- 
icality in systems with attractive interactions but there 
are some differences. The fluid and crystal states are dis- 
tinguished by their spatial periodicity and fluidity, in ad- 
dition to their density. Therefore, there is another tran- 
sition^) corresponding to the two properties in the su- 
percritical region. When periodic order is not established 
even after fluidity is lost, there must be an intermediate 
phase, which is considered the glass phase [l(| . The tran- 
sition from fluid to glass is considered to be dynamical 
transition corresponding to ergodicity breaking jlll [l2| ■ 

The phase behavior of a polydisperse HS system still 
remains under discussion. Bartlett and Warren studied 
polydisperse systems using a density functional theory 
(DFT) and claimed that the thermodynamic function 
does not have a singularity at the point of equal concen- 
tration and that the first-order transition line is extended 
to high-density region to surround the crystal phase fl3j . 
Furthermore, Fasolo et al. pointed out the importance 
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of fractionation; segregation into multiple crystals. Each 
crystal has different mean radius and relative dispersion 
is small inside it. By considering the free energies of 
mixed states, a lot of coexisting phases appears and the 
phase diagram becomes very complicated [14J| . Although 
such fractionated states is reasonable in thermal equilib- 
rium state, these phases are not observed in experiments 
or numerical simulations. One reason is that the system 
is glassy in the regime where fractionation is predicted. 
Therefore diffusion of particles over long distance, which 
is necessary for segregation, is suppressed. 

In the present paper, we don't treat fractionation 
but consider long surviving homogeneous state including 
both of equilibrium and nonequilibrium ones. Especially, 
we discuss the relation between the fluid-crystal transi- 
tion for small dispcrsity and the fluid-glass transition for 
large dispcrsity. We perform nonequilibrium molecular- 
dynamics (MD) simulations, which is not only useful to 
study nonequilibrium dynamics but also gives clue to 
reveal equilibrium property. We study the fluid-crystal 
transition in equilibrium for low dispersity by nonequilib- 
rium simulation. On this aspect, a number of numerical 
studies on polydispcrse hardcore systems has been previ- 
ously reported. But these have been highly restricted to 
two-dimensional hard-disk systems [15|, Qj| . Since two- 
dimensional systems show peculiar properties owing to 
the low dimensionality, the study of three-dimensional 
systems is necessary. On the other hand, it is difficult 
to perform simulations with sufficiently large linear di- 
mensions in three dimensions, thus finite size effect often 
makes the conclusions ambiguous. Nonequilibrium anal- 
ysis without time-consuming equilibration makes large- 
scale simulations possible. 

Let us denote the contents of the present paper. In 
the next section, the detail of the model and method of a 
numerical simulation is explained. From section III to V, 
we investigate three transitions among fluid, crystal and 
glass states to obtain the nonequilibrium phase diagram 
shown in Fig. [T] The final section is devoted for the 
concluding remarks. 



II. MODEL: HARD ELASTIC SPHERES 

We perform MD simulations of clastic spheres with 
a fixed number of particle N, temperature T and pres- 
sure P using the Nose-Hoover method [l7|, E3 and the 
Parinello- Rahman method [l9j]. The reason we did not 
employ the standard event-driven simulation of HSs (20| 
is that pressure control, which is essential in the following 
analysis of first order transition, cannot be implemented 
efficiently to this method. Since hard-sphere system is 
widely accepted as one of reference models for solid-fluid 
transitions, we estimate hard-limit of elastic modulus by 
extrapolation, which is described later. Hereafter, we use 
the units with which temperature ksT = 1, mean radius 
fi = 1, and the mass of the particle to be mi = 1. 

Polydispcrsity is introduced by a uniform distribution 



of particle radii. The strength of polydispcrsity is mea- 
sured by the standard deviation, A = \l (r^ — l) 2 , where 
Ti is the radius of particle i and 777 denotes the average 
over all particles. It is known that the quantitative prop- 
erties of a polydisperse system are well described only by 
A and that the detailed form of the distribution function 
of Ti is irrelevant when polydispcrsity is not too strong 

0- 

The interaction between contacting particles, i and j, 
is given by the Hertzian contact potential, £t)[|qi — qj| — 
( r i + r j)] 5 ^ 2 j where q, is the position of particle i. The 
interaction energy equals zero when |q^ — q.,-| > + ?*-,-. 
The system becomes a true HS when Young's modu- 
lus Eq approaches infinity. Young's modulus is set to 
Eo = 10 4 — 10 7 . Since we use finite values of Young's 
modulus, the particles are allowed to overlap to make 
the effective radius and density smaller. In order to cor- 
rect this effect, we consider the effective hardcore packing 
fraction <p\ lCl which corresponds to the density of the sys- 
tem with an infinite Young's modulus. By considering 
the equipartition of energy, the overlap length of parti- 
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cles is proportional to E Q . The packing fraction of 
the corresponding hardcore system </>h c is therefore ex- 
pected to be (/>hc = (4tt/3V) J2i — co^o^ 5 ) w i tri a 
calibration constant cq. This constant is determined to 
be Co = 0.48 by performing preliminary simulations, with 
which we confirmed that the extrapolated state equation 
exhibits good agreement with the result of event-driven 
simulation of the true HS system. This correction is used 
throughout the letter. For example, </>h c is 0.6% smaller 
than (in/W) J2i r$ for E = 10 6 . 

We adopt two types of initial particle configuration; 
FCC and random configurations. The radii of parti- 
cles are assigned randomly in accordance with the dis- 
tribution mentioned above and independently of the po- 
sitions. For random initial state, we perform simulations 
with overdamped dynamics before integrating Hamil- 
ton's equations of motion until the maximum kinetic en- 
ergy of the particles decreases below 200 to avoid the 
rapid acceleration of strongly coalesced particle pairs. 
After that, the initial velocities are randomly assigned 
by the Boltzmann distribution. 

The observed quantities discussed below are the data 
for Eq = 10 6 averaged over 4 samples with TV = 55296, 
unless otherwise stated. We confirm that our conclusions 
do not change in a larger system with N = 131072. Time 
integration is performed by the fourth-order predictor- 
corrector method using a discrete time step of At = 
0.0004 — 0.02 and typical number of integration step is 
4 x 10 6 . 

III. CRYSTAL MELTING TRANSITION 

First, we analyze the polydispcrsity dependence of the 
fluid-crystal transition and clarify the existence of the 
predicted critical point (P C ,A C ). The order parameter 
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FIG. 2: (color online) Time derivatives of the mean density 
is plotted with respect to pressure for A = 0.0 — 0.085. The 
derivative d<f>/dt is approximated by [4>{2t w ) — (j)(t w )] jt w with 
t TO /0.001 = 2 19 and 2 20 . 



corresponding to this criticality is the density gap 5<j) be- 
tween the bistable phases at the melting point. This 
is calculated by a two-step simulation. As a first step, 
we determine the melting pressure P m (A) for a given 
A(< A c ) from the nonequilibrium analysis discussed 
later. After that, we observe the packing fractions of the 
bistable states, fluid (A,P TO (A)) and so r ld (A, P m (A)), 
individually by performing equilibrium simulations with 
both fluid (random) and crystal(FCC) initial conditions. 
The packing fraction of the fluid/crystal at the melt- 
ing point gives the lower/upper bound of the coexisting 
phase at a fixed (j> condition (see Fig. [I) and its width is 
(A,P m (A)) - 0fl uid (A,P m (A)). 

To determine the melting point, we observe the 
nonequilibrium relaxation from the mixed initial state 
[2TI I22I] ; a half of the cubic space is occupied by the 
crystal and the remaining part is occupied by the ran- 
dom packing (fluid) . Thus the two regions are separated 
by a flat interface, which is perpendicular to the (100)- 
direction of the FCC structure at time t = 0. As t in- 
creases, the interface moves so that the fraction of the 
phase with lower free energy increases. The melting pres- 
sure P m can be determined as the point where the sign 
of d(j)/dt in the steady state changes, since positive and 
negative values of d(f>/dt indicate that crystallization and 
melting occur at the interface, respectively. This method 
requires a relatively short-time simulation compared to 
the equilibrium method and enables us to treat larger sys- 
tems and reduce the finite-size effect. Figure [2] shows the 
pressure dependence of dcf)/dt and obtained P m (A) gives 
the phase boundary between fluid and crystal phases in 
the inset of Fig. [1] 

By this steady interface motion, we can compare equi- 
librium stability of FCC and fluid states. There remains a 
possibility, however, that there can be more stable state, 
such as other types of crystal structure. But we expect 
that it is not the case for polydispersity, A < 0.088. 

Performing additional equilibrium simulations at these 
P m (A), we eventually obtain A-dependence of the 5(f> 
shown in Fig. [3] This order parameter approaches zero 




FIG. 3: (color online) Polydispersity dependence of the den- 
sity gap between the fluid and crystal states on the melting 
line. The results for different Young's moduli are shown to- 
gether but very little difference is observed. The solid curve 
denotes 8<j> = 0.45 x (0.088 - A) 0,73 , (inset) Log-log plot of 
the relation between the density gap and the crystalline order 
parameter on the melting line. 

as 8(j) oc (A c - A)^ 1 with A c = 0.088(2) and /3 = 0.7(2). 
The critical pressure P c = 3.0(2) and the critical pack- 
ing fraction cj) c — 0.576(4) are also obtained. We also 
calcu late the FCC order param eter of the crystal state; 
to = cos[Kpcc ' — where Kpcc is the fun- 

damental reciprocal vector of the FCC crystal. The in- 
set of Fig. [3] indicates a power-law; m oc d(j)(A)^ m ^ oc 
(A c - A) 0m with p m = 0.04(1). The range of observed 
value of to is, however, too narrow to conclude the exis- 
tence of the power-law. 



IV. GLASS TRANSITION 

We next investigate the transition between the fluid 
and glass phases by scanning P at fixed A(> A c ). Sig- 
nificant change is observed around a certain threshold 
P 9 (A); the mobility of particles markedly decreases ap- 
proaching P g , which denotes the glass transition point. 

Figure 0] sho ws the P-depcndenc e of the diffusion con- 
stant, D(t w ) = \<\i(2t w ) — qi(t w )\ 2 /t w , with waiting time 
t w under the random initial condition. Here we make 
the time interval to measure the displacement equiva- 
lent with t w so that only single time scale is introduced. 
While the D(t w ) converges to a certain equilibrium value 
by increasing t w for P < P g w 3.0, the relaxation is so 
slow for P > P g that equilibrium state cannot be ob- 
tained for used values of t w . Instead, we remarks on 
the aging property, i.e., the persistent waiting-time de- 
pendence; D(t w ) continues to decrease with t w , roughly 
in a power law, above P g , This indicates that as the 
relaxation proceeds, the system becomes trapped in an 
increasingly stable metastable state and the dynamics 
becomes slower. 

As clearly observed in Fig. [H the D{t w ) vs P curves 
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FIG. 4: (color online) Pressure dependence of the diffusion 
constant for fixed polydispersity. The data for various waiting 
times are plotted together to show the aging behavior. The 
initial state is random packing with packing fraction 0.50. 
(inset) Pressure dependence of the packing fraction for fixed 
polydispersity. For each A, we show the data at three times, 
i/0.001 = 2 20 ,2 21 , and 2 22 to show good convergence. 



hardly depend on A for A < 0.12. Thus -P g (A) is also in- 
dependent of A and equals to 3.0, similar to the value of 
P c . In general, the effect of polydispersity is small except 
in the crystal phase. In addition, almost the same behav- 
ior is observed even in the subcritical region (A < A c ), 
as a supersaturation phenomena fl2| . Any sign of crystal 
nucleation is not observed for A > 0.60. It is known that 
polydispersity drastically reduces the nucleation rate of 
the crystal Q. In the inset of Fig.[H 0hc is plotted with 
respect to P. The packing fraction also has little de- 
pendence on A (slightly increases with A) both in the 
fluid and glass phases. Therefore, the glass transition 
density <j> g (A) = (j)(P g (A)) also has little dependence 
on A and <f> g ps 0.57 ps <p c . The extrapolated value, 
9 (A — > 0) ps 0.57, agrees with the value, 4>d ps 0.58, pre- 
dicted by mode-coupling theory [llj or mean field theory, 
which corresponds to the appearance of the exponentially 
many metastable states in the fluid (23j . 

Above the threshold pressure P g , 4>h c gradually ap- 
proaches the random close packing (RCP) fraction 
0rcp(A) (24[, which equals 0.635 for the monodisperse 
(A = 0) system and increases very slowly with A [25| . 



V. CRYSTAL-AMORPHOUS TRANSITION 

Finally, we consider the transition between the crys- 
tal and glass states, which is driven by sweeping A at 
fixed P(> P c )- In Fig- El we plot the A-dependcncc 
of 0hc, estimated by simulations under both crystal and 
random (glass) initial conditions. Both the crystal and 
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FIG. 5: (color online) Polydispersity dependence of the pack- 
ing fraction at fixed pressure. The data for FCC and ran- 
dom initial configurations are plotted together, which coin- 
cide for large A. In each case, we show the data at three 
times, t/0.001 = 2 16 ,2 17 , and 2 18 . 



glass states are (meta)stable for a long time in the high- 
pressure region since there is too little free volume for the 
local structures to reconfigurate. The packing fraction of 
the crystal is larger than that of the glass for small A but 
they becomes indistinguishable above a certain threshold 
A am (P). This threshold appears to be universal, i.e., it 
hardly depends on P. In addition, its value is very sim- 
ilar to A c ps 0.088. The density difference continuously 
decreases to zero as A approaches A am . Similar behavior 
is observed for the crystalline order parameter under the 
FCC initial condition. 



VI. CONCLUSIONS 

In summary, we investigated the dynamical transitions 
of a polydispcrsc clastic sphere system by MD simula- 
tions by remarking on noncquilibrium states including 
metastable states. The obtained transition lines with 
respect to the packing fraction and polydispersity are 
summarized in Fig. Q] It was confirmed that the first- 
order transition between the fluid and crystal phases ter- 
minates at the critical point (<^> c , A c ) and that the other 
two phase boundaries begin from the critical point to sur- 
round the glass phase. The glass state has intermediate 
properties between those of the fluids and crystal states; 
it exhibits temporal freezing but does not have periodic 
order. The fluid-glass and crystal-glass boundaries can be 
drawn in a surprisingly simple way and are expressed as 
P ~ P c (or <fi ~ <fi c ) and A ~ A c , respectively. The glass 
transition line passes through the critical point, which is 
reasonable because the continuous breakdown of the crys- 
tal requires marginal fluidity at the critical point. While 
softening of the interaction potential will not make essen- 
tial change in the phase behavior, that of a system with 
attractive interactions is an interesting open problem. 

The transition between the fluid and glass states is 
not considered to be an equilibrium transition but a dy- 
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namical one since the static quantities do not exhibit 
any singular behavior and the transition line is elongated 
into the crystal phase in equilibrium. This continuous 
relationship from the supercritical region to the super- 
saturating monodisperse system suggests the equivalence 
of the dynamical glass transitions in the monodisperse 
and poly disperse systems. Let us note that the criti- 
cal packing fraction is close to the random loose parking 
(RLP) fraction, <^rxp ~ 0.56, which is considered to be 
the minimum packing fraction required to maintain the 
internal stress for highly frictional particles [26}. This 
coincidence seems natural considering that RLP gives a 
criterion related to the excluded volume effect. The free 
volume of particles is very small above (/>rlp and diffusion 
is highly suppressed. 

Let us consider the meaning of the boundary between 
the crystal and glass states. The transitions at this 
boundary are continuous in terms of density, in contrast 
to those of fluid-crystal boundary in the subcritical re- 
gion. It is natural that the first-order transition line and 



a continuous transition line should meet at the multi- 
critical point. But this is conflict with the prediction of 
first order transition by Bartlett and Warren [l3j. Our 
nonequilibrium analysis cannot eliminate the possibility 
that the crystal phase is metastable below A am , which is 
estimated by nonequilibrium simulations, and first-order 
transition occurs at A < A am . We wonder, however, 
whether a mcan-ficld-likc approach in the DFT scheme 
can treat the state around the terminal point, where the 
fluid exhibits singular behavior in dynamics and the pe- 
riodicity of the crystal is damaged. Our numerical result 
suggests the possibility that criticality remains. 
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